classdef TrussElement < Element
    % TRUSSELEMENT 桁架单元类
    % 用于三维桁架结构，只承受轴向力
    
    properties (Constant)
        ElementType = 'Truss'
        NumNodes = 2           % 桁架单元只有2个节点
        Dimension = 3          % 支持三维空间
        DofPerNode = 3         % 每个节点3个平动自由度
    end
    
    properties (SetAccess = protected)
        Length                 % 单元长度
        Orientation           % 单元方向向量
    end
    
    methods
        function obj = TrussElement(connectivity, material, integrationRule)
            % 构造函数
            if nargin == 0
                connectivity = [1, 2];
            end

            if nargin <= 2
                integrationRule = GaussIntegrationRule(1); % 桁架单元只需1点积分
            end
            
            obj@Element(connectivity, material, integrationRule);      
        end
        
        function Ke = computeStiffnessMatrix(obj)
            % 计算桁架单元刚度矩阵 (三维)
            obj.validateConfiguration();
            
            % 获取单元属性
            L = obj.computeLength();              % 单元长度
            D = obj.MaterialModel.getConstitutiveMatrix();

            % 局部坐标系下的刚度矩阵 (1D)
            k_local =  D / L * [1, -1; -1, 1];
            
            % 转换到全局坐标系
            T = obj.computeTransformationMatrix();
            Ke = T' * k_local * T;
        end
        
        function Me = computeMassMatrix(obj, massType)
            % 计算桁架单元质量矩阵
            if nargin < 2
                massType = 'consistent';
            end
            
            obj.validateConfiguration();
            
            density = obj.MaterialModel.Density;
            A = obj.MaterialModel.Area;
            L = obj.computeLength();
            totalMass = density * A * L;
            
            if strcmpi(massType, 'lumped')
                % 集中质量矩阵 - 质量平均分配到节点
                Me = (totalMass / 2) * eye(6); % 每个节点1/2质量
            else
                % 一致质量矩阵
                Me = obj.computeConsistentMassMatrix(density, A, L);
            end
        end
        
        function Fe = computeLoadVector(obj, loadCase)
            % 计算桁架单元载荷向量
            obj.validateConfiguration();
            
            nDofs = obj.getTotalDofs();
            Fe = zeros(nDofs, 1);
            
            switch lower(loadCase.Type)
                case 'bodyforce'
                    % 考虑自重等体力
                    Fe = obj.computeBodyForceVector(loadCase.Value);
                otherwise
                    error('TrussElement: Unsupported load type: %s', loadCase.Type);
            end
        end
        
        function stress = computeStress(obj, nodalDisplacements)
            % 计算桁架单元轴向应力
            obj.validateConfiguration();
            
            % 获取材料属性
            E = obj.MaterialModel.YoungsModulus;
            L = obj.computeLength();
            
            % 计算轴向应变
            axialStrain = obj.computeAxialStrain(nodalDisplacements);
            
            % 计算轴向应力 (胡克定律)
            axialStress = E * axialStrain;
            
            stress = struct();
            stress.axial = axialStress;
            stress.vonMises = abs(axialStress); % 桁架单元von Mises应力就是轴向应力绝对值
            stress.strain = axialStrain;
        end
        
        function strain = computeStrain(obj, nodalDisplacements)
            % 计算桁架单元轴向应变
            obj.validateConfiguration();
            strain = obj.computeAxialStrain(nodalDisplacements);
        end
        
        function N = evaluateShapeFunctions(obj, naturalCoords)
            % 计算桁架单元形函数 (一维自然坐标 ξ ∈ [-1,1])
            xi = naturalCoords;
            N = 0.5 * [1 - xi, 1 + xi];
        end
        
        function dN_dxi = evaluateShapeGradients(obj, naturalCoords)
            % 计算形函数导数 (对自然坐标ξ)
            dN_dxi = 0.5 * [-1, 1]; % 常数
        end
        
        function J = computeJacobian(obj)
            % 计算雅可比矩阵 (桁架单元是常数)
            L = obj.computeLength();
            J = L / 2; % 一维雅可比
        end
        
    end
    
    methods (Access = protected)
        function [B, detJ] = assembleStrainDisplacementMatrix(obj, naturalCoords)
            % 组装桁架单元的应变-位移矩阵B
            % 对于桁架单元，B矩阵是常数
            
            L = obj.computeLength();
            dirCos = obj.computeOrientation();
            
            % 转换到全局坐标系
            B = zeros(1, 6); % 1个应变分量，6个自由度
            B(1, 1:3) = -dirCos / L;   % 节点1的贡献
            B(1, 4:6) = dirCos / L;    % 节点2的贡献
            
            detJ = L / 2;
        end
        
        function L = computeLength(obj)
            % 计算单元长度
            if isempty(obj.NodeCoordinates)
                L = 0;
                return;
            end
            
            node1 = obj.NodeCoordinates(1, :);
            node2 = obj.NodeCoordinates(2, :);
            L = norm(node2 - node1);
        end
        
        function dirCos = computeOrientation(obj)
            % 计算单元方向余弦向量
            if isempty(obj.NodeCoordinates)
                dirCos = [1, 0, 0]; % 默认方向
                return;
            end
            
            node1 = obj.NodeCoordinates(1, :);
            node2 = obj.NodeCoordinates(2, :);
            dirVec = node2 - node1;
            dirCos = dirVec / norm(dirVec);
        end
        
        function T = computeTransformationMatrix(obj)
            % 计算局部到全局坐标的转换矩阵
            dirCos = obj.computeOrientation();
            
            % 对于三维桁架，转换矩阵为 [l, m, n] 的扩展形式
            T = zeros(2, 6);
            T(1, 1:3) = dirCos;   % 节点1的方向余弦
            T(2, 4:6) = dirCos;   % 节点2的方向余弦
        end
        
        function axialStrain = computeAxialStrain(obj, nodalDisplacements)
            % 计算轴向应变
            L = obj.computeLength();
            dirCos = obj.computeOrientation();
            
            % 提取节点位移
            u1 = nodalDisplacements(1:3)'; % 节点1位移
            u2 = nodalDisplacements(4:6)'; % 节点2位移
            
            % 计算轴向位移差在单元方向上的投影
            du = (u2 - u1) * dirCos';
            
            % 轴向应变 = 轴向位移差 / 原长
            axialStrain = du / L;
        end
        
        function Me = computeConsistentMassMatrix(obj, density, A, L)
            % 计算一致质量矩阵
            totalMass = density * A * L;
            
            % 局部坐标系下的一致质量矩阵
            Me_local = (totalMass / 6) * [2, 1; 1, 2];
            
            % 转换到全局坐标系
            T = obj.computeTransformationMatrix();
            Me = T' * Me_local * T;
        end
        
        function validateConfiguration(obj)
            % 验证桁架单元配置
            if isempty(obj.NodeCoordinates)
                error('Node coordinates must be set before computation');
            end
            if isempty(obj.MaterialModel)
                error('Material model must be set before computation');
            end
        end
    end
    
end